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The numerical simulation of strongly first-order phase transitions has remained a notoriously difficult problem 
even for classical systems due to the exponentially suppressed (thermal) equilibration in the vicinity of such a 
transition. In the absence of efficient update techniques, a common approach to improve equilibration in Monte 
Carlo simulations is to broaden the sampled statistical ensemble beyond the bimodal distribution of the canon- 
ical ensemble. Here we show how a recently developed feedback algorithm can systematically optimize such 
broad-histogram ensembles and significantly speed up equilibration in comparison with other extended ensem- 
ble techniques such as flat-histogram, multicanonical or Wang-Landau sampling. As a prototypical example of 
a strong first-order transition we simulate the two-dimensional Potts model with up to Q — 250 different states 
on large systems. The optimized histogram develops a distinct multipeak structure, thereby resolving entropic 
barriers and their associated phase transitions in the phase coexistence region such as droplet nucleation and 
annihilation or droplet-strip transitions for systems with periodic boundary conditions. We characterize the effi- 
ciency of the optimized histogram sampling by measuring round-trip times t(N, Q) across the phase transition 
for samples of size TV spins. While we find power-law scaling of r vs. TV for small Q < 50 and TV < 40 2 , 
we observe a crossover to exponential scaling for larger Q. These results demonstrate that despite the ensem- 
ble optimization broad-histogram simulations cannot fully eliminate the supercritical slowing down at strongly 
first-order transitions. 

PACS numbers: 02.70.Rr, 64.60.De, 75.10.Hk 



I. INTRODUCTION 

Competing phases or interactions in many-particle systems 
can give rise to complex free energy landscapes that exhibit 
multiple local minima and maxima. Thermal equilibration in 
these systems slows down exceedingly due to the (exponen- 
tial) suppression of tunneling across free energy barriers. Ex- 
amples of such slowly equilibrating systems can be found in 
various settings ranging from spin glasses to folded proteins. 

Numerical approaches to simulate these systems in ther- 
mal equilibrium suffer from the same slowing down when 
based on variational techniques or conventional Monte Carlo 
sampling. To overcome this bottleneck various alternative 
sampling approaches have been developed in recent years. 
Most of these approaches, which include multicanonical sam- 
pling HIEI, broad-histogram sampling [3|, parallel temper- 
ing 1UJ6), multiple Gaussian modified ensemble Q, and 
Wang-Landau sampling EJHJ, can be broadly grouped as ex- 
tended ensemble techniques. Their common goal is to sample 
a statistical ensemble that allows to significantly broaden the 
range of sampled energies beyond the comparatively narrow 
distribution of the canonical ensemble. 

The Wang-Landau algorithm tries to bring the idea of a 
broad sampling to an extreme by sampling a flat histogram 
in energy space. However, it was soon realized that sam- 
pling a uniform energy distribution is not necessarily the opti- 
mal way to improve equilibration and reduce autocorrelation 
times lfT0l - [T2l . Instead it turns out that in order to (consid- 
erably) speed up equilibration and minimize autocorrelation 
times one should sample a non-uniform energy distribution 
that allocates more statistical weight to the bottleneck(s) of 



the simulation which typically coincide with the free energy 
barriers lfT3l . These so-called optimized ensembles are tai- 
lored to a given physical system and directly reflect the un- 
derlying free energy landscape. One can systematically ob- 
tain the optimized statistical ensemble from an initial broad- 
histogram distribution by applying a feedback algorithm lfL3ll 
that reallocates statistical weight based on measurements of 
the (local) diffusivity of the random walk which the system 
performs in energy space during the simulation. This ensem- 
ble optimization has been applied in a broad variety of phys- 
ical systems suffering from long thermal equilibration times 
in the absence of efficient non-local updates including folded 
proteins lfT4HT7l . frustrated magnetic systems |fl8l |T9ll , and 
dense liquids l20l . The technique has further been general- 
ized to optimize the grid of temperature points used in paral- 
lel tempering simulations [21 J, has been used in combination 
with cluster udpates [22] and has been adopted for the simu- 
lation of quantum systems 112311241 . 

In this manuscript, we apply and analyze the ensemble opti- 
mization technique in the context of a strong first-order phase 
transition where the characteristic double-well structure of the 
free energy provides a generic situation for entropically sup- 
pressed equilibration. In particular, we consider the thermal 
phase transition of the Q-state Potts model in the limit of large 
Q, with our calculations being performed for up to Q = 250 
states. We find that the optimized ensemble aims to overcome 
the entropic barrier(s) of this transition by allocating most of 
the statistical weight in the energy range that corresponds to 
phase coexistence, e.g. the suppressed energy region of the 
characteristic bimodal distribution of the canonical ensemble. 
Remarkably, a multi-peak structure evolves in the optimized 
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phase space 

Figure 1: (color online) Sketch of the free energy landscape in phase 
space for a slowly equilibrating system. 



histogram that clearly resolves various intermediate transi- 
tions between metastable states, such as droplet formation and 
droplet-strip transitions. 

The remainder of the manuscript is structured as follows: 
We will first provide a brief review of the ensemble optimiza- 
tion technique in Section [II] In the subsequent Section III 



we 



turn to the large-Q Potts model and discuss the multiple dis- 
tinct features of the optimized broad-histogram distribution. 
We conclude our analysis by measuring the performance of 



the optimized ensemble technique in Section IV 



II. OPTIMIZED ENSEMBLES 

We start our discussion of the ensemble optimization tech- 
nique by first offering a broader view on Monte Carlo sam- 
pling and statistical ensembles. We then briefly review the 
derivation of the optimized ensembles and a related feedback 
algorithm. 



A. Monte Carlo sampling and statistical ensembles 

Speaking in broader terms one might take the perspective 
that the idea underlying Monte Carlo sampling is to map a 
random walk in some high-dimensional space of configura- 
tions {Cj} 



ci 



C2 



Ci+l 



onto a random walk in a lower dimensional space, such as 
energy space (which is a one-dimensional space) 



E( Cl ) -> E(c 2 ) 



E{c t ) -> E(a +1 ) 



and to define a statistical ensemble in this latter low- 
dimensional space which then determines the transition 
probabilities between configurations in the original high- 
dimensional space. In particular, the statistical ensemble as- 
signs a statistical weight w(c) to a configuration c solely on 
the basis of the respective energy E(c) of that configuration 

w(c) = w(E(c)) . 



The most commonly used statistical ensemble, of course, is 
the canonical ensemble, where the statistical weights are de- 
fined as 

w{c) = exp(-PE(c)) . 

In order to simulate a reversible Markov process in configura- 
tion space one then defines transition probabilities from con- 
figuration c to c such that detailed balance is fulfilled. Com- 
mon choices for these transition probabilities are Metropolis 
weights 

^Metropolis C -> c ) = mm 1) 7~T = mm 1, , 

V w(c)J V w{E{c))J 



or heat -bath weights 



Pheat-bath(c -> £) 



w{c) 



Km) 



w(c) + w(c) w(E(c)) + w(E(c)) ' 



While these choices of transition probabilities indeed ensure 
that the random walk in configuration space is Markovian, 
it should be noted that the projected random walk in energy 
space is not Markovian. This becomes clear when consider- 
ing that multiple configurations may have the same energy E, 
whereas the distribution of energies that can be reached by a 
single update may be completely different for each of these 
configurations. Thus, there is additional information encoded 
in configuration space which is not captured by E, and it is 
this 'memory' which makes the projected random walk in en- 
ergy space non-Markovian. 



B. Non-uniform diffusivity and optimized histograms 

The random walk in energy space has another distinct fea- 
ture: the (local) diffusivity of this random walk, which for a 
given energy level measures the ability of the random walker 
to move to other energy levels, is not uniform in energy space. 
In fact, it is exactly this modulation of the local diffusivity 
which reflects the roughness of the underlying energy land- 
scape. A suppressed diffusivity signals a 'bottleneck' of the 
simulation and is typically associated with a phase transition 
or other entropic barrier. 

This modulation of the local diffusivity thus differentiates 
the various energy regimes for a given system, and in this 
light it becomes clear that one shortcoming of flat-histogram 
techniques is that they use a uniform distribution of statisti- 
cal weight across these inherently different energy regimes. 
In contrast the optimized ensemble method allocates statisti- 
cal weight based on measurements of the local diffusivity and 
shifts additional statistical weight towards the bottleneck(s) of 
the simulation, e.g. those energy regimes with a suppressed 
local diffusivity. As a result the so-optimized random walk 
in energy space will sample a non-uniform histogram, spend 
more time in energy regimes with low diffusivity, and thereby 
do its best to suppress the bottlenecks associated with the un- 
derlying free energy landscape. 

In more technical terms, we consider a random walk in 
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some energy range [E_, E + ] between two extremal energies 
E- and E + . In this paper we sample the entire energy range, 
so E- and E + are, respectively, the lowest and highest ener- 
gies that our model has. The random walkers in energy space 
will drift between these two extremal energies and we can 
think of the overall random walk as being composed of two 
opposite steady-state 'currents' between these two extremal 
energies. These two currents exactly compensate one another, 
as the system remains in equilibrium, and are independent of 
energy. We can express these currents as 



j = D(E)H(E) 



4f_ 

dE 



(1) 



where D(E) is the local diffusivity in energy space, H(E) is 
the sampled energy histogram and f(E) defines the orienta- 
tion of the current by measuring for a given energy the frac- 
tion of random walkers which have last visited one of the two 
extremal energies, say, the lower extremal energy E_. This 
latter fraction can be measured by recording two histograms, 
H + (E) and H_(E), where, for each Monte Carlo step, one 
increments the histogram with label '+' or '-' depending on 
which extremal energy the random walker has visited last. 
The two histograms H + (E) and if_ (E) thus sum up to the to- 
tal histogram H{E) = H+{E) + H_(E). The fraction f(E) 
is then given by f(E) = H_ (E) /H(E). 

In order to speed up equilibration one wants to maximize 
the current ([TJ between the two extremal energies. Varying 
the histogram H(E) this can be achieved ITT31 by sampling a 
non-uniform distribution 



H° V \E) 



1 



y/D(E) 



(2) 



that is inversely proportional to the square root of the local dif- 
fusivity D(E) and thus reallocates statistical weight to those 
energy levels with suppressed diffusivity. 



C. The feedback algorithm 



In order to sample the optimized ensemble (|2]i we apply 
the feedback algorithm outlined in Ref . [ 1 3 1 . We start from 
an initial broad-histogram ensemble with statistical weights 
w(E) which we obtain from a few iterations of the Wang- 
Landau algorithm or by extrapolating results from smaller 
system sizes. Running a (short) simulation for this initial en- 
semble we record the two histograms H + (E) and H_ (E) in- 
troduced above which in turn allow us to calculate the local 
diffusivity as 



D(E) cx H{E) 



'dE 



(3) 



We then refine the statistical weights by feeding back this local 
diffusivity and define new optimized weights as 



Wopt(E) =w(E)a 



1 1 df 
H(E) dE 



(4) 



Subsequent simulations are performed for this new set of sta- 
tistical weights. To further improve and eventually converge 
the statistical weights for the optimized ensemble we repeat 
the feedback procedure several times. Note that in order to 
ensure convergence the number of Monte Carlo steps between 
subsequent feedback iterations needs to be increased; we typi- 
cally double the number of Monte Carlo steps for consecutive 
runs. 



D. Improving the first feedback step 

There is a certain trade-off in performing the early feed- 
back steps in the algorithm outlined above: On the one hand, 
an early feedback after only a small number of Monte Carlo 
sweeps appears advantageous as it may quickly give dra- 
matically improved statistical weights and thereby speed up 
all subsequent simulations. On the other hand, the quality 
of the feedback is rather sensitive to noisy input data, es- 
pecially in calculating the numerical derivative df /dE used 
in the feedback. To minimize this trade-off one thus needs 
a way to quickly estimate this latter derivative in the pres- 
ence of (substantial) noise. Conventional approaches such as 
finite-difference formulas, however, turn out to be exquisitely 
sensitive to the noise in the recorded histograms H + (E) 
and H_{E). In particular, the measured fraction f(E) — 
H-(E)/ H(E) is a monotonically decaying function only 
when the simulation is in equilibrium in the simulated sta- 
tistical ensemble which for a suboptimal choice, such as the 
flat-histogram ensemble, may require rather long Monte Carlo 
runs. 

We have therefore developed a scheme that allows for the 
estimation of the derivative in the presence of significant 
noise. The idea is to analyze the measured fraction f(E) in 
Fourier space, truncate the high-frequency terms which can be 
associated with noise, and then determine the derivative using 
the low-frequency terms only. In doing so, we make use of the 
fact that for a continuous Fourier transformation 



i 
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f(E)dE, 



the derivative of the original function f(E) 
1 



d E f(E) 



J- 



luj ■ e 



f(uj) doj 



(5) 



can be easily calculated in Fourier space and then transformed 
back. 

In implementing this idea one needs to work around several 
obstacles. First, in order to avoid irrelevant boundary terms, 
the function to be analyzed using the Fourier transformation 
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should be periodic. We therefore concatenate f(E) with its 
reflection. Secondly, the above relation strictly holds only for 
the continuous Fourier transformation. As the energy levels 
of the Potts model and, in general, the energy bins of a broad 
histogram simulation are discrete, we need to work with a dis- 
crete Fourier transformation. To overcome errors introduced 
by this, we make use of the following iterative scheme which 
refines the calculated derivative by iteratively reducing the de- 
viation between the integral of the approximated derivative 
and the original function: 



a) internal energy 





Sf 1 = B E f(E) 
Sf = Sf 1 + ~d E [ f(E) 



SfdE 



Sf +L = 5f + d E [f(E)- / SfdE 



Here, Be denotes the approximate derivative using the 
Fourier-based scheme above. The scheme is iterated until the 
norm of the correction term falls below a certain threshold. 



IH. THE LARGE-Q POTTS MODEL 

The two-dimensional Q-state Potts model is well 
known E5ll to undergo a thermal phase transition which 
turns from continuous for small Q < 4 to weakly first-order 
for Q > 4 and eventually becomes a strong first-order 
transition for Q 3> 5. We will turn to this latter case of a 
strong first-order transition in systems with up to Q = 250 
different Potts states to explore the extent to which the 
optimized ensemble algorithm can achieve equilibration at 
such a transition. 

The Hamiltonian of the Q-state Potts model is given in 
terms of spins a t which take discrete values a t £ {1, . . . , Q} 
as 



(6) 



Here the sum runs over all pairs of nearest neighbors on a 
square lattice, and the Kronecker (5-function tests whether two 
Potts spins have the same values. We have run simulations for 
two sample geometries: a "toroidal geometry", i.e. a square 
lattice with periodic boundary conditions, and a "cube geom- 
etry" by forming a cube with square lattices on each of its 6 
faces. 

We will start our discussion by briefly mentioning both ex- 
act and numerical results for thermodynamic properties of the 
Potts model in this large, but finite Q limit. We will then turn 
to the energy regime associated with phase coexistence at this 
first-order phase transition and examine the various intermedi- 
ate, metastable states such as droplets or strips which occur in 
this regime. In particular, we will discuss a distinct multi-peak 
structure which emerges in the optimized histogram distribu- 
tion and show how these features can be linked to transitions 
between the various metastable states. 
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Figure 2: (color online) Thermodynamic properties of the Q-state 
Potts model in the large, but finite Q limit: a) energy, b) specific 
heat, c) free energy, and d) entropy. The temperature axis is rescaled 
by the transition temperature T* = 1/ln (1 + y/Q). Data shown is 
for system size L = 22 with periodic boundary conditions. 
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A. Thermodynamic properties 

The thermal phase transition in the Q-state Potts model 
from a disordered phase at high temperatures to an ordered 
phase at low temperatures occurs at a transition temperature 
T* which for the infinite system is found |25 1 to be 



ln(l + y/Q) 



(7) 



This phase transition is accompanied by sharp features in vari- 
ous thermodynamic properties such as the energy, the specific 
heat, the free energy and the entropy. Since the optimized en- 
semble algorithm allows to directly calculate the density of 
states g(E), we can readily compute all of these thermody- 
namic variables. Our results are summarized in Fig.|2]for sim- 
ulations with Q = 10, 50, 250 states, where we have rescaled 
the temperature axis by the transition temperature T* in the 
thermodynamic limit. As expected, the features associated 
with the phase transition sharpen as we increase the number 
of Potts states Q for a system of fixed size L. For instance, 
the discontinuous jump in the energy grows with increasing Q 
and U(T) approaches a step function in the limit of Q — » oo. 
It is this broadening energy regime within the discontinuous 
jump of the energy that is associated with phase coexistence 
and the occurrence of intermediate, metastable states as we 
will discuss in detail below. 



B. Phase coexistence and metastable states 

The distinct characteristic of a first-order phase transition is 
a free energy profile that passes through a double well shape 
as one drives the transition with some external parameter such 
as temperature. At the transition temperature the two minima 
of the free energy are exactly equal leading to coexistence of 
the two phases. For the high-Q Potts model at hand there is a 
considerable amount of latent heat associated with this transi- 
tion, i.e. the internal energies of the two phases in proximity 
to this phase transition vastly differ as shown in Fig. [2^). As 
the system goes from one phase to the other this latent heat 
is not released (or absorbed) in a single step, but the system 
undergoes a sequence of phase transitions between various 
metastable states which are not minima of the free energy in 
thermal equilibrium, but correspond to states with intermedi- 
ate internal energies. One such metastable state is a droplet of 
one phase inside the other phase. Since the free energy density 
of the two phases becomes arbitrarily close in the vicinity of 
the transition, the free energy cost of forming a droplet is due 
to the surface of the droplet, and not to its volume. It is thus 
entropically favorable to nucleate and grow a single droplet of 
a shape that minimizes its surface free energy. This droplet 
condensation transition has recently been studied in detail for 
a variety of physical systems using both numerical ifTUl l26l - 
28 1 and analytical [29-33 1 approaches. 

For a torus geometry, e.g. a system with periodic bound- 
ary conditions, this droplet will subsequently expand as the 
total energy is changed until it percolates and it becomes en- 



tropically more favorable to form a strip wrapping around one 
of the boundaries. As this strip further grows the role of the 
two phases will eventually be reversed and the system will un- 
dergo a second sequence going from a strip to a droplet and 
eventually annihilate the remaining droplet to complete the 
phase transition from one phase to the other. This transition 
was first discussed for the Ising model in an external magnetic 
field by Leung and Zia [34 1 and studied in detail by Neuhaus 
and Hager [ 10 1 using multicanonical Monte Carlo sampling. 
The droplet nucleation and droplet-strip transitions were also 
observed for the Potts model with Q = 10 and system sizes of 
up to 1024 x 1024 spins using a microcanonical approach l35l . 

We show representative snapshots of spin configurations re- 
flecting these metastable states in Fig. [3] All snapshots have 
been taken taken from our numerical simulations of the 250- 
state Potts model. 



C. Droplet nucleation and droplet-strip transition 

For the toroidal system we can thus distinguish four inter- 
mediate transitions taking place "within" a first-order transi- 
tion: The nucleation of a dominant droplet (which might oc- 
cur via the condensation of multiple small droplets), a droplet 
to strip transition and two more processes where the roles of 
the two phases are reversed. These intermediate transitions 
occur at energies that are within the discontinuous jump of 
the internal energy U{T) plotted in Fig. [3] a) and are therefore 
not equilibrium energies at any temperature. For the canon- 
ical ensemble the states at these intermediate transitions are 
strongly suppressed, with its characteristic bimodal distribu- 
tion of sampled energies as shown in the bottom panel of 

Fig- a 



1. Multi-peak structure 

In sharp contrast to the canonical ensemble the feedback al- 
gorithm of Sec. reallocates significant statistical weight to 
the energy range located within the double-peak structure of 
the canonical distribution corresponding to the discontinuous 
jump of the internal energy. Strikingly, we find the emergence 
of a distinct multi-peak structure in this energy range as shown 
in the top panel of Fig. [4] The emergent peaks resolve pre- 
cisely the four intermediate transitions discussed above. We 
come to this identification, as given in Fig. [4] by i) compar- 
ing the energies of typical configuration snapshots as shown 
in Fig. [3] with the locations of these peaks, ii) estimating the 
transition energies of the droplet-strip transitions as discussed 



in Section nn~C2 and iii) calculating order parameters for the 



droplet-strip transitions as detailed in Section III C 3 The re- 
distribution of statistical weight in this multi-peak structure 
also reveals that these transitions between metastable states 
are of different severity. With the histogram peaks corre- 
sponding to the droplet-strip transitions being much more pro- 
nounced than those corresponding to droplet nucleations we 
can conclude that the entropic barriers associated with the 
former transitions are significantly larger than those associ- 
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a) E/2N = -0.20 b) E/2N = -0.35 c) E/2N = -0.54 d) E/2N = -0.87 




Figure 3: (color online) Snapshots of spin configurations in the phase coexistence region of the 250-state Potts model, a) Formation of several 
small ordered droplets within the disordered phase, b) A dominant droplet is formed, c) The droplet percolates to a strip wrapping around one 
boundary, d) A single disordered droplet remains in the ordered phase. All data shown is for system size L = 20. 



ated with the latter transitions. Another observation regarding 
this emerging multi-peak structure is that the histogram distri- 
bution is not perfectly symmetric with respect to the ordered 
/ disordered phases. For instance, the difference of the two 
smaller peaks reflects that droplet formation in the disordered 
phase is associated with a larger entropic barrier than droplet 
formation in the ordered phase. 

These characteristic features of the multi-peak structure fur- 
ther evolve as we vary the strength of the underlying first- 
order transition by increasing the number of Potts states Q 
or the system size L as shown in Figs. [5] and [6] respectively. 
With increasing the number of Potts states Q we find the 
droplet-strip transitions to attract considerably more statistical 
weight than the droplet formation transitions. In particular, 



droplet strip 
transition 




-1 -0.8 -0.6 -0.4 -0.2 



energy E I 2N 

Figure 4: (color online) Histograms of the optimized ensemble (top 
panel) and the canonical ensemble at the transition temperature T* 
(bottom panel) for the 250-state Potts model with L = 22 and 
toroidal geometry. In contrast to the bimodal distribution of the 
canonical ensemble the histogram of the optimized ensemble reveals 
a distinct four-peak structure reflecting the transitions between the 
various metastable states in the phase coexistence region. 
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Figure 5: (color online) Optimized histograms for increasing num- 
ber of Potts states Q and fixed system size L = 22 with toroidal 
geometry. While the histogram peaks in the emerging multi-peak 
structure seem to diverge with increasing Q for intermediate ener- 
gies and associated with the droplet-strip transitions, the histogram 
peaks associated with droplet formations appear to converge. The 
plotted energy range corresponds to the coexistence region defined 
in Eq. (jTTJ. 

the histogram peaks associated with the droplet-strip transi- 
tions seem to diverge with increasing Q, while the histogram 
peaks associated with the droplet formation transitions appear 
to converge to a finite height while sharpening with increasing 
Q, see Fig. [3] Similarly, we find that increasing the system 
size L increases the peaks associated with the droplet-strip 
transitions more strongly than those associated with the for- 
mation of a droplet, as shown in Fig. [6] 



2. Location of droplet-strip transitions 

We can estimate the location of the intermediate droplet- 
strip transitions more quantitatively by estimating the inter- 
face length of the droplet / strip on either side of the tran- 
sition. Making such an estimate for the droplet, however, 
requires knowledge of its rough shape. The latter depends 



7 




tions can be approximated as 



-0.6 -0.4 
rescaled energy E ' 



Figure 6: (color online) Optimized histogram for increasing sys- 
tem size L and fixed number of Potts states Q = 50. Similar to 
Fig.[5]the peaks associated with the droplet-strip transitions prolifer- 
ate more strongly with increasing system size than those associated 
with droplet formation. The plotted energy range corresponds to the 
coexistence region defined in Eq. \\\\ . 



on the anisotropy of the surface tension and to some extent 
the geometry of the system. For the Potts model, the surface 
tension and consequently the equilibrium droplet shape have 
been calculated analytically for a droplet of fixed size embed- 
ded in an infinite volume [36 37 1; for finite systems, Billoire 
et al. |38] have estimated the surface tension based on multi- 
canonical simulations of the canonical probability density for 
mixed-phase states. In the limit of large Q, which we consider 
here, the surface tension a is found to become isotropic and 
the droplet shape is expected to be roughly circular. The lo- 
cation of the droplet-strip transition can then be estimated by 
identifying the radius of the droplet R for which the surface 
free energy of the droplet, i*dropiet = IttuR, becomes equal to 
the surface free energy of the strip, -F s trip = 2erL. The transi- 
tion is therefore expected to take place at R — L/ir. At this 
transition point, the droplet occupies area ttR 2 = L 2 /ir and 
thus a fraction 1/tt of the total area of the sample. We can 
estimate the total energy of a given domain pattern as a sum 
of the contributions from the domains and the interfaces. For 
the example of an ordered droplet of radius R in a disordered 
background we have 



E=(L 2 - 7ri? 2 )e dis (T) + irR 2 e old (T) + 2irRe a 



(8) 



This configuration, by local stability of the curved interface, 
is at a temperature T that is below the transition temperature 
T* by an amount proportional to the curvature 1/R of the 
interface. At that temperature, the energy densities of the two 
phases are e or d(T) and edi s (T), while e a is the excess energy 
per unit length in the interface. 

Simplifying this expression by keeping only the contribu- 
tions that are proportional to the areas of the domains, we ob- 
tain that the transition energies of the two droplet-strip transi- 



E. 



droplet-strip, 1 



E, 



droplet-strip, 2 



1 

7T 

= -1- 



(9) 
(10) 



where these energies are given relative to the size of the co- 
existence region, e.g. the width of the jump of the internal 
energy plotted in Fig. [3^). In the following, we rescale en- 
ergies such that the energy of the ordered phase E or( i(T*) is 
mapped to — 1 at the transition temperature T* and the energy 
of the disordered phase Edi s (T*) becomes 0: 



E* 



E - E md (T*) 
E dis (T*) - E olA (T* 



(11) 



To rescale our numerical results we use the exactly-known en- 
ergy densities of the two phases at the transition temperature 
T* in the thermodynamic limit. 

We have indicated the so-estimated locations of the two 
droplet-strip transitions by the vertical bars in Figs. [5] and [6] 
Indeed, the respective histogram peaks associated with these 
transitions seem to converge to these locations in the limit of 
large Q and L. In more quantitative terms, the energy of the 
interface moves E* at the transitions to a higher energy by 
an amount proportional to 1/L. The small shift in T due to 
the curvature of the interface of the droplet also moves E* 
at the transition by amount proportional (ignoring logs) to 
l/(Ly/Q) at large Q; this latter effect is of the same sign at 
the lower-energy droplet-strip transition, where the system is 
slightly "superheated". The trends with increasing Q and L 
at this lower transition can be clearly seen in Figs. [5] and [6] 
respectively. At the higher-energy droplet-strip transition the 
two 1/L finite-size effects are of opposite signs, so the peak 
in the histograms stays closer to E* = —1/tt. 

At the energies of these droplet-strip transitions, the two 
configurations (droplet and strip) have the same entropy. 
However for the system to make a transition between these 
two configurations, it must increase the amount of interface 
by an amount proportional to L. The entropy deficit per unit 
length in the interface is proportional to logQ at large Q, 
so the entropy barriers at these transitions are proportional to 
L log Q. 



3. Order parameter for the droplet-strip transition 

An alternative approach to locate and further describe the 
droplet-strip transition is to define an order parameter for this 
transition. In doing so we follow an idea of Neuhaus and 
Hager [10| and measure the existence of a strip by measur- 
ing the dimensions Li,L 2 of the minimal bounding box for 
the largest droplet in the system 



dis/ord = § r L max ( L d 



dis/ord T-dis/ord 



(12) 



where the index "dis/ord" distinguishes whether the phase in 
the droplet corresponds to the disordered or ordered phase, re- 
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Figure 7: (color online) Optimized histogram and susceptibilites for 
the droplet percolation order parameters. The lower susceptibility 
peak has been rescaled to match the height of the peak in the op- 
timized histogram. Data shown is for a 250-state Potts model and 
system size L = 24. 



spectively. When the droplet percolates and a strip is formed, 
one dimension of the bounding box becomes equal to the sys- 
tem size L and the order parameter jumps from to 1. Fur- 
thermore, we can associate a susceptibility with this order pa- 
rameter, 



xo = (O 2 ) - (O) 2 = (O) - (O) 2 e [0,1/4] 



(13) 



which we expect to proliferate at the droplet-strip transition. 
Comparing the divergence of this susceptibility with the re- 
spective intermediate peaks forming in the optimized his- 
togram, we find that they coincide not only in location, but 
also their respective shapes as shown in Fig. [7] The latter illus- 
trates that the entropic barriers at this intermediate transition 
which the optimized ensemble overcomes by shifting statisti- 
cal weight towards this transition arise solely from percolating 
a droplet into a strip. 
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Figure 8: (color online) Anisotropy of the largest droplet of ordered 
phase as denned in Eq. {14} in the phase coexistence region of a 
250-state Potts model (L = 18). A spherical droplet corresponds to 
anisotropy 1 as indicated by the dashed line. At the droplet-strip tran- 
sitions indicated by the vertical bars the droplet anisotropy quickly 
changes. The deviation from unity above the upper droplet-strip tran- 
sition indicates an anisotropy due to finite-size effects. 



For energies above the upper droplet-strip transition we find 
that the droplet anisotropy deviates from unity which indicates 
that the ordered droplet in an otherwise disordered phase ex- 
hibits an non-spherical rather than a spherical shape. Since 
the analytical calculation for a droplet embedded in an infinite 
system predicts a roughly spherical shape 071 . the observed 
anisotropy must be rooted in the finite size of our system. 
Interestingly, we also seem to observe a signature indicating 
the droplet formation with the droplet anisotropy suddenly in- 
creasing from a rj 1.25 to a w 1.5 in the respective energy 
region. 



D. Simulations on the cube surface 



4. Droplet anisotropy 

Finally, we return to the droplet anisotropy and estimate 
corrections to the circular shape induced by the finite size of 
our system. To this end we monitor the anisotropy of the 
droplet by measuring the ratio of the dimensions L\,Li of 
the minimal bounding box for the largest (ordered) droplet in 
the system 

«= m ^44- (14) 

In this notation a spherical droplet corresponds to a = 1. 

In Fig. [8] we plot this anisotropy for the largest ordered 
droplet measured for energies in the phase coexistence re- 
gion. At the droplet-strip transitions - indicated by the verti- 
cal bars in the figure - the droplet anisotropy undergoes rapid 
changes as the droplet percolates and deforms into a strip. 



Since the droplet-strip transition causes the main bottle- 
neck for simulations in a toroidal geometry, one could ask 
whether simulations on other surface topologies, in particular 
on simply-connected surfaces, do not suffer from entropic bar- 
riers originating from such shape transitions. In order to ex- 
plore this idea we have simulated the Potts model in a "sphere 
topology" by considering the surface of a cube. We assem- 
ble such a cube surface with L sites on each edge and a total 
number of sites N = 6L 2 — 12L + 8 such that the corners 
have coordination number z = 3, while all other sites have 
coordination number z = 4. 

Fig. [9] shows results for the optimized histograms for the 
Potts model on such a cube surface. In panel a), histograms 
for fixed Q = 250 and L in the range L = 6 ... 16 are shown, 
while in panel b), the system size is fixed at L = 10 and 
simulations are shown for various Q = 10 . . . 250. As op- 
posed to the simulations on the torus, the peaks correspond- 
ing to droplet nucleation/evaporation are now dominating the 
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Figure 9: (color online) Optimized histogram for the Potts model on 
the cube surface for (a) a fixed number of Potts states Q = 250 and 
various system sizes L, (b) a fixed system size L — 10 and different 
choices of the number of Potts states Q. In contrast to the simulations 
on the torus (cf. Fig.[6|, the dominant peaks are related to droplet nu- 
cleation/annihilation transitions, while the emerging, smaller peaks 
reveal transitions between states with droplets occupying an increas- 
ing number of corners of the cube. The plotted energy range corre- 
sponds to the coexistence region defined in Eq. {TTJ. 



histogram for these parameters. 

However, for the largest systems with Q = 250 and L > 14 
(N > 1016 spins) spins, four smaller peaks emerge. Exam- 
ining snapshots of the system in the associated energy ranges, 
we find that droplets nucleate around corners of the cube and 
the transitions mark a change in the number of occupied cor- 
ners. A similar observation has been made previously for 
multicanonical simulations of the Ising model in an external 
field Q0). 

There are multiple reasons why droplets nucleate at the cor- 
ners of the cube. Naively, one might think that this is solely 
due to the lower coordination number of the corner sites, 
which provides a small energetic advantage (of order one in 
the system size) to place a droplet on a corner. Closer in- 
spection of the surface free energy of a droplet enclosing a 
corner, however, reveals that similar to the droplet-strip tran- 
sition there are entropic barriers which scale with the system 
size L. To see this consider a droplet of area A. If A is small 
enough (relative to L 2 ), this droplet may sit on one of the faces 
of the cube so that it encloses no corners of the cube. Then the 
surface it sits on is flat, so the surface free energy is minimized 
by a circular droplet with radius R so that A = irR 2 . Alter- 
natively, the drop may enclose one corner. Putting the corner 
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Figure 10: (color online) Comparison between optimized histograms 
for the cube surface and the torus at Q — 50 for various system sizes 
L. The plotted energy range corresponds to the coexistence region 
defined in Eq. {TTJ. 



at the center of the droplet, the droplet can be a quarter-circle 
on each of the 3 adjacent faces. These quarter-circles each 
have radius R with A = (37ri? 2 )/4. The net result is that the 
perimeter of the droplet (which sets the surface free energy) 
is smaller by a factor of v3/2 when the droplet is centered 
on a corner as compared to when it does not contain a corner. 
The difference is of order y/A, so for large enough cubes will 
dominate over the order one effect mentioned above. 

While the simulation bottlenecks / entropic barriers asso- 
ciated with the corner transitions are significantly suppressed 
in comparison with the droplet-strip transitions of the toroidal 



geometry, which is illustrated Fig. 10 we have not succeeded 
in suppressing all entropic barriers of order L by going from 
the torus to the cube surface. As a consequence, we expect 
the same asymptotic performance of the optimized ensemble 
simulations for both geometries. 



IV. SAMPLING EFFICIENCY AND ROUND-TRIP TIMES 

We finally address the numerical efficiency of sampling 
the optimized statistical ensemble for the Q-state Potts model 
and, more generically, for a strong first-order phase transi- 
tion. In order to quantify this performance we follow ear- 
lier work lfTTI - [T3l and measure the characteristic time scale of 
the random walk in energy space, e.g. the round-trip time to 
traverse the full energy range [£L, E + ], and its scaling with 
system size L and number of Potts states Q. Our results are 



summarized in Figs. 11 and 12 



For systems undergoing continuous transitions it was 
shown 1 1 1 1 that the flat-histogram ensemble sampled in the 
Wang-Landau method generally suffers from a power-law 
slow down 



7"flat-histogram (X N 2 ■ N Z 



(15) 
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Figure 11: (color online) Scaling of the round-trip time between the 
lowest and highest energy states with system size L for the optimized 
ensemble simulations using heat bath dynamics and the toroidal ge- 
ometry. While the scaling can be fitted to polynomial behavior for 
small Q < 35, the scaling crosses over to exponential behavior for 
larger Q > 50. 



Figure 12: (color online) Scaling of the round- trip time between the 
lowest and highest energy states with the number of Potts states Q for 
the optimized ensemble simulations using both Metropolis and heat 
bath dynamics in the toroidal geometry. The dashed lines represent 
cubic regressions to the data. 



which is reminiscent of the well-known critical slowing down 
in the canonical ensemble. The additional exponent z depends 
on the system at hand, and was measured to be z ps 0.4 for 
the Ising model and z ps 0.9 for the fully frustrated Ising 
model ifTTl . In contrast, the optimized ensemble method does 
not suffer from such a 'critical slowing down' at continuous 
transitions and produces round-trip times that scale almost 
perfectly with system size 



^optimized— ensemble (X -^V ' (In iV) 



(16) 



up to a logarithmic correction |13|. The improved scaling of 
the optimized ensemble can thus considerably speed up the 
sampling efficiency of a broad-histogram simulation, with im- 
provements of two orders of magnitude reported already for 
intermediate system sizes [ 1 3 1 . 

Exploring the efficiency of simulations at a first-order phase 
transition, previous work has found an exponential divergence 
of the round-trip time with L in the Ising model in an ex- 
ternal magnetic field [ 10 1, which is referred to as supercrit- 
ical slowing down. This divergence occurs in the presence 
of shape transitions such as the droplet-strip transition dis- 
cussed above and the quantitative behavior of the scaling can 
indeed be related to the surface tension of droplets. In a recent 
study Neuhaus et al.. [7 1 demonstrated that the multiple Gaus- 
sian modified ensemble (MGME) method suffers only from a 
residual supercritical slowing down when applied to the Potts 
model (with Q = 20 and Q — 256 states in their simulations). 

For the optimized ensembles the scaling of the round-trip 
times is shown in Figs. [TT|and[T2|as a function of both system 
size L and number of Potts states Q, respectively. Similar to 
previous results for continuous transitions, we find a dramatic 
improvement of the optimized ensemble when compared to 
the flat-histogram ensemble. However, even for the optimized 
ensemble we do not recover the almost perfect scaling ( |T6] > 
observed for continuous transitions when the severity of the 



first-order transition proliferates, e.g. by increasing the sys- 
tem size L or the number of Potts states Q. In general, we 
find that the scaling is independent of the transition dynamics, 
e.g. whether we choose Metropolis or the heat-bath transition 
rates. 

We first look at the scaling with system size L when fix- 
ing the number of Potts states Q in the range 8 < Q < 150. 



As shown in Fig. 1 1 we find that, for the range of L studied, 
the round-trip times appear to scale polynomially r ~ N 2+z 
in system size L for Q < 35. We estimate the effective ex- 
ponents to be z ps 0.3 for Q = 8, z ps 0.31 for Q = 10, 
z fa 0.38 for Q = 20, and z ps 0.48 for Q = 35. We 
note that in this regime the first-order transition is still rel- 
atively weak and we do not (yet) observe the characteristic 
multi-peak structure in the optimized histogram. However, 
when further increasing number of Potts states Q we find that 
the scaling crosses over to the expected exponential behavior 
r ~ exp(aL) for Q > 50 due to the entropy barriers at the 
droplet-strip transitions, and this becomes apparent already on 
rather small system sizes. This, of course, bears witness of the 
strong first-order nature of the phase transition in this regime, 
which also becomes evident in a noticeable multi-peak struc- 
ture of the optimized histogram even for considerably small 
system sizes L > 14. Our results however do not exclude a 
crossover to supercritical scaling also for weaker transitions 
at some larger system size. Indeed, results in Ref. 35 indicate 
that also for Q = 10, an extended stripe phase emerges for 
system sizes L > 512, which may well lead to the same effect 
we observe for large Q at much smaller L. 



V. CONCLUSIONS 

We presented an application of the optimized ensemble 
method to improve equilibration of broad-histogram simula- 
tions of the first-order transition in the large-Q Potts model. 
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The optimized histogram develops a characteristic multipeak 
structure which indicates that the system releases latent heat in 
a sequence of phase transitions. The intermediate metastable 
states exhibit droplets of the coexisting phases in varying 
shapes. For a toroidal system geometry the dominant bottle- 
neck in the simulation is the entropic barrier associated with a 
droplet-strip transition. 

We find that the ensemble optimization is capable to only 
partially overcome this bottleneck by shuffling statistical 
weight towards the entropic barriers. It still exhibits the 
same asymptotic supercritical slowing down for strongly first- 
order transitions previously reported for flat-histogram simu- 
lations [10]. It thus remains an open question whether simula- 
tion schemes for strongly first-order phase transitions can be 



further improved to overcome this asymptotic slowing down. 
Possibly, statistical ensembles defined in multiple system vari- 
ables could further improve extended ensemble simulations, 
or one might turn to modified update technique, which, for 
instance, attempt to specifically update the boundaries of a 
droplet. 
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